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Abstract  - The  stability  of  two-dimensional 
Finite-Difference  Time-Domain  subgridding  schemes 
was  numerically  examined.  Both  the  same-time-step  and 
the  multiple-time-step  schemes  were  considered.  Results 
show  that  the  multiple-time-step  subgridding  scheme  is 
late-time  unstable  due  to  larger-than-unity  eigenvalues. 
As  to  the  same -time-step  subgridding  schemes,  stability 
is  related  to  the  treatment  of  comer  regions. 

Keywords  - FDTD  method,  subgridding,  stability. 

I.  INTRODUCTION 

Since  the  Finite-Difference  Time-Domain  (FDTD) 
method  was  introduced  [1],[2],  a major  challenge  is 
modeling  locally  fine  structures  and/or  field  singularities. 
To  model  fine  features  or  resolve  rapid  field  variations, 
small  cells  are  required  and  the  overall  computational 
cost  can  be  prohibitive  by  using  the  traditional  FDTD 
method. 

One  attractive  solution  is  subgridding,  i.e.,  fine 
meshes  are  employed  wherever  needed  and  the  rest  of 
the  computational  domain  is  discretized  by  coarse 
meshes.  The  early  version  of  subgridding  applies  the 
FDTD  method  twice  in  different  regions,  where  the 
results  in  coarse  region  are  used  as  boundary  conditions 
for  fine  regions  [3].  Due  to  the  lack  of  timely  feedback 
from  fine  regions,  this  approach  only  yields  better  results 
in  fine  regions.  To  improve  the  accuracy,  fine  regions 
are  embedded  into  coarse  regions  and  the  timely 
feedback  is  provided  by  field  coupling  on  boundaries  [4]. 
In  this  scheme,  two  methods  are  possible  to 
accommodate  the  stability  requirement  imposed  by  fine 
region  cell  size.  One  is  to  use  the  same  time  step  (STS) 
size,  which  is  dictated  by  the  fine  region  cell  size,  in  both 
regions.  The  other  is  to  use  different  time  step  sizes 
determined  by  coarse  and  fine  region  cell  sizes 
respectively.  Thus  one  time  step  in  coarse  region 
corresponds  to  multiple  time  steps  (MTS)  in  fine  region. 
Accordingly,  temporal  interpolation  is  required  for  field 
coupling  on  the  boundaries. 

Despite  the  various  accuracy-enhancement  procedures 
[5]-[8],  FDTD  subgridding  methods  are  notorious  for 


their  late-time  instability,  which  typically  happens  after 
many  iterations.  To  address  this  issue,  spatial  reciprocity 
on  coarse-fine  region  boundaries  has  proposed  as  a 
necessary  condition  to  achieve  late-time  stability  [9]. 
However,  its  effectiveness  remains  unclear.  In  this  paper, 
we  provide  numerical  examinations  of  its  sufficiency  for 
both  the  STS  and  the  MTS  subgridding  schemes. 

The  dominant  eigenvalues  of  the  system  amplification 
matrices  are  examined  as  a key  measure  of  stability. 
Meanwhile,  corresponding  numerical  simulations  are 
carried  out  for  sufficient  time  steps  to  check  the  actual 
late-time  behavior.  For  simplicity,  2D  TEz  wave  was 
examined.  The  subgridding  refinement  ratio  is  3:1.  The 
odd  ratio  results  in  collocated  coarse-fine  region 
boundary  values  in  space  and  in  time.  For  other 
subgridding  ratios,  the  analysis  follows  the  same 
procedure.  We  note  that  the  main  purpose  of  this  paper  is 
to  provide  first-hand  numerical  examinations. 
Theoretical  explanations  are  given  tentatively  based  on 
known  theories. 

This  paper  is  organized  as  follows.  Section  II  presents 
the  methodology  and  describes  the  numerical  tests. 
Section  III  verifies  the  approach  by  the  regular  FDTD 
scheme.  Sections  IV  and  V show  the  results  of  the  STS 
and  the  MTS  schemes,  respectively.  Finally,  conclusion 
remarks  are  drawn  in  section  VI. 

n.  METHODOLOGY 

The  stability  of  FDTD  subgridding  schemes  was  first 
studied  by  the  dominant  eigenvalues  of  the  system 
amplification  matrices  [2].  In  the  discretized  Maxwell’s 
equations 


[A]  denotes  the  system  amplification  matrix.  To  verify 
the  late-time  behavior,  simulations  were  carried  out 
inside  Perfect  Electric  Conducting  (PEC)  boxes  of 
various  sizes.  A magnetic  point  source  near  the 
lower-left  comer  was  excited  and  the  electric  field  at 
another  location  was  monitored  at  each  time  step.  The 
source  is  a differentiated  Gaussian  pulse  and  its  peak 
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magnitude  is  less  than  one.  Therefore,  the  observed 
electrical  field  magnitude  should  not  be  larger  than  one 
in  a stable  simulation.  Whenever  the  above  condition  is 
violated,  we  claim  that  a simulation  is  unstable  and  that 
time  step  is  recorded  as  the  “critical  point”.  Otherwise, 
the  most  recent  observations  were  checked  after  ten 
million  time  steps.  If  they  do  not  show  any  sign  of 
growth,  a scheme  is  claimed  as  numerically  stable. 

The  fine  region  is  placed  at  the  center  of  the  com- 
putational domain,  which  is  denoted  by  m *n,  where  m 
is  the  distance  from  the  fine  region  boundary  to  the  outer 
boundary  of  the  entire  computational  domain  and  n the 
size  of  the  fine  region,  both  are  in  terms  of  the  coarse 
region  cell  size. 

We  note  that  the  maximum  possible  value  at  the 
observation  is  different  with  respect  to  different  com- 
putational domain  sizes.  The  “critical  point”  only 
indicates  instability  and  its  value  should  not  be  in- 
terpreted in  the  same  way  for  different  computational 
domains.  Due  of  limited  resources,  we  only  considered 
small  computational  domains.  However,  stability  is 
usually  characterized  asymptotically  by  assuming  an 
infinitely  large  one.  To  examine  the  differences  and 
verify  our  approach,  we  first  study  the  regular  FDTD 
method. 

IH.  THE  REGULAR  FDTD  METHOD 

Figure  1 shows  the  dominant  eigenvalues  of  the 
system  amplification  matrix  of  the  regular  FDTD 
method  as  a function  of  the  Courant-Friedrichs-Levy 
(CFL)  number,  which  is  defined  by  A/V2vp  / Ah  • Here,  vp 

is  the  phase  velocity,  Ah  is  the  cell  size  and  At  is  the 
time-step  size.  The  computational  domain  size  is  de- 
noted by  m *n  in  terms  of  cells,  where  m represents  the 
length  and  n represents  the  width. 

It  is  interesting  to  notice  that  simulations  of  very  small 
computational  domains  can  be  stable  even  with  a larger 
than  one  CFL  number.  For  a 3 x 3 PEC  box,  the  largest 
stable  CFL  number  reaches  to  1.15  (the  CFL  numbers 
are  resolved  to  0.01  in  this  article).  As  the  computational 
domain  size  increases,  the  largest  stable  CFL  number 
drops  rapidly  and  becomes  one  when  the  computational 
domain  size  is  equal  to  or  larger  than  12  x 12.  These 
results  were  confirmed  by  simulations. 

The  above  observations  can  be  understood  by  the 
asymptotic  nature  of  theoretical  stability  analysis,  which 
only  considers  the  worst  case  in  the  Von  Neumman 
analysis  [2].  It  is  well  known  that  the  time  step  size  of  2D 
FDTD  simulations  is  restricted  by 


where  kx  and  ky  are  associated  with  the  discrete  Fourier 
modes  [10]  and  the  worst  case  is  kxh  = kyh  = (2n+l)n.  In 
an  infinitely  large  computational  domain,  there  are 
infinite  Fourier  modes  including  the  above  worst  case. 
When  the  computational  domain  size  is  limited,  the 
discrete  Fourier  modes  that  can  possibly  exist  are  limited 
and  they  may  not  include  the  worst  case.  The  stability 
condition  determined  directly  by  the  dominant 
eigenvalues  of  small  computational  domains  can  be 
different  from  that  of  the  Von  Neumman  analysis. 
Therefore,  multiple  tests  with  increasing  computational 
domain  sizes  are  suggested  to  estimate  the  asymptotic 
behavior  of  numerical  stability. 


Fig,  1.  Dominant  eigenvalues  of  the  system 

amplification  matrix  of  the  regular  2D  FDTD 
method  according  to  different  CFL  numbers  and 
computational  domain  sizes. 


IV.  THE  STS  SUBGRIDDING  SCHEME 

In  general,  we  can  write  the  STS  scheme  as 


^ «+l/2  ^ 

(A 

ncc 

-n+U2 

\U  J 

kAfc 

A ff) 

'UH 

where  U 0+1/2  contains  all  coarse  region  unknowns  and 

u”*1 i 2 contains  all  fine  region  unknowns,  both  at  the 
current  time  step.  ACc  represents  how  the  coarse  region 
uses  the  previous  coarse  region  values,  ACf  represents 
how  the  coarse  region  uses  the  previous  fine  region 
values,  AFC  represents  how  the  fine  region  uses  the 
previous  coarse  region  values  and  AFF  represents  how 
the  fine  region  uses  the  previous  fine  region  values. 

Figures  2 and  3 depict  a regular  coarse-fine  region 
boundary  and  a comer  region  respectively.  In  order  to 
calculate  the  fine  region  boundary  electric  field,  e.g.  eyIi 
we  need  the  coarse  region  magnetic  field,  e.g.  h’zI.  To 
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recover  the  missing  values  due  to  boundary  truncation, 
the  following  linear  interpolation  is  employed  [8], 

K,  =\h„+\h»  • (4) 

<s> 

K,=Hzl.  (6) 


Instead  of  being  half-cell  away  from  the  boundary 
electric  field,  the  recovered  magnetic  field  is  1.5-cell 
away.  To  keep  at  least  first-order  accuracy  in  space,  we 
employed  the  unbalanced  differencing  scheme  to  update 
the  boundary  electric  field,  e.g., 


e 


n 

yi 


At 

+ 2A he 


kr112 


(7) 


where  z//?  is  the  fine  region  cell  size.  When  calculating 
the  coarse  region  boundary  magnetic  field,  e.g.  Hz2 , we 
need  the  electric  field  along  the  coarse-fine  region 
boundary,  e.g.  Ey2.  They  are  obtained  by  enforcing 
spatial  reciprocity  [9],  i.e., 


+ ey4j+  e 


y 3 


- (8) 


Note  that  the  coarse  region  boundary  magnetic  field  is 
calculated  by  central  differencing  scheme  in  the  layout 
shown  in  Fig.  2.  Another  choice  is  to  extend  the  fine 
region  1/3  coarse  cell  (or  one  fine  cell)  to  the  left,  i.e., 
overlapping  the  coarse  and  fine  regions  by  one  fine  cell. 
In  that  case,  the  coarse  region  boundary  magnetic  field  is 
calculated  by  unbalanced  differencing  and  the  fine 
region  boundary  electric  field  is  calculated  by  central 
differencing.  Since  unbalanced  differencing  scheme  is 
first-order  accurate,  choosing  fine  region  boundary 
electric  field  to  be  calculated  by  unbalanced  differencing 
scheme  is  apparently  more  accurate. 

As  we  shall  see,  the  real  difficulty  is  to  treat  comers, 
e.g.,  h’zl  and  in  h”zJ  Fig.  3,  by  imposing  spatial 
reciprocity.  If  we  write 


*;.=!».■ +5^.^=!^ +\n«. 


(9) 


the  coarse  region  boundary  electrical  field  can  be  ob- 
tained by  spatial  reciprocity  as, 

= ^jc3  = ^ y2  + 

This  apparently  introduces  errors  when  calculating  Hz3 
(also  HzI  and  in  Hz2).  To  improve  the  accuracy,  other 
comer  layouts  can  be  applied  by  shifting  the  fine  region 
boundary  field  components  and  violating  the  reciprocity. 
Since  spatial  reciprocity  is  our  major  concern,  we 
employ  the  layout  shown  in  Fig.  3 and  refer  to  it  as 
“reciprocal  with  comer”  treatment. 
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Fig.  2.  The  coarse-fine  region  boundary. 


Fig.  3.  The  coarse-fine  region  boundary  at  a comer. 


Alternatively,  one  may  simply  ignore  comers  and  write 

Kl=Hli,K>=Hl2.  (11) 

This  can  be  physically  interpreted  as  treating  Hzi  and  Hz2 
as  the  average  in  those  cells.  In  the  following,  we  refer  to 
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this  as  “reciprocal  without  comer”  treatment. 

A.  Reciprocal  with  Corner 

An  example  of  the  system  amplification  matrix  is 
illustrated  in  Fig.  4,  which  corresponds  to  a “3  x 3” 
subgridding  region.  The  dominant  eigenvalue  of  this 
example  is  1.420754  when  the  CFL  number  is  1.03.  As 
we  decrease  the  CFL  number,  the  dominant  eigenvalue 
drops  quickly  and  becomes  1.0  when  the  CFL  number  is 
1.01.  We  further  calculated  the  dominant  eigenvalues  of 
subgridding  schemes  with  larger  computational  domain 
sizes  (up  to  “10  x 10”).  It  was  found  that  the  results  are 
all  1 .0  when  the  CFL  number  is  1 .0. 

Although  the  dominant  eigenvalues  do  not  show  any 
sign  of  instability  for  CFL=1,  numerical  simulations 
were  unstable  for  all  computational  domains  being 
tested.  Moreover,  the  instability  does  not  happen  late. 
For  example,  the  “3  x 3"  case  became  unstable  after  360 
time  steps.  The  “9  x10"  case  is  the  most  stable  one, 
which  only  runs  stably  for  3689  time  steps.  This  example 
shows  that  enforcing  spatial  reciprocity  on  coarse-fine 
region  boundaries  is  not  a sufficient  condition  for 
stability  in  general. 


COLUMN 

Fig.  4.  An  example  of  the  system  amplification  matrix  of 


the  “reciprocal  with  comer”  treatment. 

B.  Reciprocal  without  Corner 

The  same  set  of  subgridding  layout  was  examined. 
The  system  amplification  matrices  look  similar  to  Fig.  4. 
The  dominant  eigenvalue  of  the  “3  x 3”  case  is  1.420752 
when  the  CFL  number  is  1.03.  As  we  decrease  the  CFL 
number,  the  dominant  eigenvalue  drops  quickly  and 
becomes  1.0  when  the  CFL  number  is  1.01.  For 
subgridding  schemes  with  larger  computational  domain 


sizes  (up  to  “10  xlO”)  the  dominant  eigenvalues  are  all 
1 .0  when  the  CFL  numbers  are  1 .0. 

Contrary  to  the  previous  case,  numerical  simulations 
were  stable  with  ten  million  time  steps  for  all 
computational  domains  being  tested. 

V.  THE  MTS  SUBGRIDDING  SCHEME 

The  major  difference  between  the  STS  and  the  MTS 
subgridding  schemes  is  the  timing  procedure  involved  in 
the  latter  case.  Figure  5 illustrates  the  case  where  one 
coarse-region  time-step  corresponds  to  three  fine-region 
time-steps.  When  calculating  en  and  en+1/s  in  fine  region, 
we  need  H”~1/6  and  }f+I/6  in  coarse  region  for  spatial 
interpolation.  Also,  when  calculating  ]~r+J/2  by  imposing 
reciprocity,  we  need  en  in  fine  region.  These  two 
requirements  make  accurate  temporal  interpolation,  i.e., 
obtaining  lf1/6  and  H”+1/6  by  IT'112  and  ff+1/29  difficult 
to  implement  if  special  reciprocity  has  to  be  enforced. 
One  solution  is  to  use  temporal  extrapolation,  i.e.,  ob- 
taining  Ifl/6  and  Hn+,/6  by  /T,/2  and  Hn'3/\  etc. 
Alternatively,  one  may  also  use  H”',/2  for  H”'l/6  and 
H” ' 1,6  [9],  which  is  frequently  referred  to  as  Zeroth 
Order  Hold  (ZOH)  in  digital  signal  processing  [11]. 
ZOH  inevitably  introduces  spectrum  distortion  to  input 
signals.  However,  we  chose  ZOH  because  temporal 
extrapolation  is  more  prone  to  instability  in  practice. 

The  implementation  of  the  above  timing  procedure 
results  in  the  following  set  of  equations, 

un~l/6  = AFFun~U2  + AFCUn'112 , (12) 

un+U6  = A2FFun-y 2 + (Aff  + 1)- AfcU"-U2,  (13) 
w”+1/2  = A3ffu"-U2  + (A2FF  + Aff  + 1)-  AfcU"-U2  ,(14) 

Un+V2  = AccUn~m  + Af-pH"-1'6 

= {Acc  + ACF-AFC)tJn-V2  (15) 
+Acf  • Affu" 

where  all  symbols  have  the  same  meaning  as  in  equation 
(3).  Accordingly,  the  system  amplification  matrix  is 
written  as 

^ _ Acc  + Acf  * Afc  Acf  • Aff  ^ 
_( ' AFF  + Aff  +\)-  Afc  Aff 

Figure  6 illustrates  the  amplification  matrix  of  the  “3 
x 3”  case.  Figure  7 shows  the  dominant  eigenvalues  vs. 
the  CFL  numbers  for  different  sizes  of  computational 
domain  and  subgridding  region.  As  we  see,  the  dominant 
eigenvalues  are  always  larger  than  1.0  regardless  of  the 
CFL  number.  As  the  size  of  computational  domain 
increases,  the  dominant  eigenvalue  for  a given  CFL 
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number  decreases  but  never  reaches  1.0.  Numerical 
simulations  further  verified  the  above  results,  where  the 
recorded  “critical  points”  are  shown  in  Fig.  8.  These 
results  correspond  to  the  dominant  eigenvalues  in  Fig.  7 
well. 

m-i  E" 

* I L ' > I I I 

itn~h  = 

Fig.  5.  Timing  of  the  MTS  subgridding  scheme,  where 
superscripts  denote  the  time-step. 


Fig.  6.  An  example  of  the  system  amplification  matrix  of 
the  MTS  subgridding  scheme. 

The  MTS  scheme  is  more  efficient  than  the  STS 
scheme.  As  demonstrated  in  Fig.  7,  the  dominant 
eigenvalues  decrease  with  either  an  increasing  com- 
putational domain  size  or  a decreasing  CFL  number. 
Thus  in  practice,  one  may  delay  the  late-time  instability 
by  either  decreasing  the  CFL  number  or  by  increasing 
the  computational  domain  size. 


Fig.  7.  Dominant  eigenvalues  of  the  MTS  subgridding 
scheme. 


Fig.  8.  “Critical  points”  of  the  MTS  subgridding 
scheme  with  different  computational  domain 
sizes  and  different  CFL  numbers. 

VI.  CONCLUSIONS 

We  numerically  examined  the  stability  of  both  the 
STS  and  the  MTS  FDTD  subgridding  schemes.  For  the 
STS  schemes,  it  was  shown  that  enforcing  spatial 
reciprocity  does  not  guarantee  stability  in  general, 
especially  when  comers  must  be  handled  in  a reciprocal 
manner.  As  to  the  MTS  subgridding  schemes,  the  system 
is  unstable  due  to  eigenvalues  that  are  out  of  the  unit 
circle.  Some  practical  considerations  were  also  given 
with  regard  to  the  use  of  FDTD  subgridding  schemes. 
Future  work  involves  developing  stable  FDTD 
subgridding  schemes. 
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